EDA + Preprocessing

data <- read.csv('data/SeoulBikeData.csv',fileEncoding = "latin1",check.names = F)
glimpse(data)
## Rows: 8,760
## Columns: 14
## $ Date                        <chr> "01/12/2017", "01/12/2017", "01/12/2017", …
## $ `Rented Bike Count`         <int> 254, 204, 173, 107, 78, 100, 181, 460, 930…
## $ Hour                        <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, …
## $ `Temperature(°C)`           <dbl> -5.2, -5.5, -6.0, -6.2, -6.0, -6.4, -6.6, …
## $ `Humidity(%)`               <int> 37, 38, 39, 40, 36, 37, 35, 38, 37, 27, 24…
## $ `Wind speed (m/s)`          <dbl> 2.2, 0.8, 1.0, 0.9, 2.3, 1.5, 1.3, 0.9, 1.…
## $ `Visibility (10m)`          <int> 2000, 2000, 2000, 2000, 2000, 2000, 2000, …
## $ `Dew point temperature(°C)` <dbl> -17.6, -17.6, -17.7, -17.6, -18.6, -18.7, …
## $ `Solar Radiation (MJ/m2)`   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
## $ `Rainfall(mm)`              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `Snowfall (cm)`             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Seasons                     <chr> "Winter", "Winter", "Winter", "Winter", "W…
## $ Holiday                     <chr> "No Holiday", "No Holiday", "No Holiday", …
## $ `Functioning Day`           <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", …
data <- data |> janitor::clean_names()
data <- data %>%
  filter(functioning_day == "Yes") |> select(-functioning_day)
data <- data %>%
  drop_na()
data <- data %>%
  mutate(holiday = ifelse(holiday == "No Holiday", 0, 1))
data <- data %>%
  mutate(time =  case_when(
    hour >=4 & hour <11 ~ "sang",
    hour >=11 & hour <13 ~ "trua",
    hour >=13 & hour <18 ~ "chieu",
    .default = "toi"
  ))
data$seasons <- factor(data$seasons)
data$holiday <- factor(data$holiday)
data$hour <- factor(data$hour)

Visualization

data |> select(-c(date,hour,seasons,holiday,time)) |>  gather(variable, value) |> group_by(variable)|> summarise(mean = round(mean(value),4),sd = sd(value),median =median(value))
## # A tibble: 9 × 4
##   variable                     mean      sd  median
##   <chr>                       <dbl>   <dbl>   <dbl>
## 1 dew_point_temperature_c    3.94    13.2      4.7 
## 2 humidity_percent          58.1     20.5     57   
## 3 rainfall_mm                0.149    1.13     0   
## 4 rented_bike_count        729.     642.     542   
## 5 snowfall_cm                0.0777   0.444    0   
## 6 solar_radiation_mj_m2      0.568    0.868    0.01
## 7 temperature_c             12.8     12.1     13.5 
## 8 visibility_10m          1434.     609.    1690   
## 9 wind_speed_m_s             1.73     1.03     1.5
data_m <-data |> select(-c(date)) |>  gather(variable, value,-c(rented_bike_count,hour,seasons,holiday,time))
ggplot(data_m,aes(x = value, y = after_stat(density))) + 
    facet_wrap(~variable, scales = "free") + 
    geom_histogram(aes(y = after_stat(density)), color = "white", fill="lightblue", alpha=1, position="identity", add_density = TRUE) + geom_density() + theme_minimal()

Tổng quan về dữ liệu:

  • Dew point temperature: Dữ liệu phân bố từ -30°C tới 25°C. Trong đó, phần lớn giá trị nằm trong khoảng nhiệt độ từ 0°C đến 25°C.

  • Humidity percent: Biểu đồ histogram cho thấy phân bố độ ẩm phần trăm từ 0% đến 100%. Các giá trị độ ẩm thường xuất hiện nhất rơi vào khoảng từ 40% đến 80%.

  • Rainfall: Lượng mưa chủ yếu tập trung ở gần giá trị 0. Dựa vào biểu đồ, có thể thấy các giá trị mưa lớn là rất hiếm trong tập dữ liệu.

  • Snowfall: Dữ liệu phân bố từ 0cm tới 9cm. Tuy nhiên, phần lớn dữ liệu tập trung ở giá trị 0cm.

  • Solar radiation: Phân bố bức xạ mặt trời tập trung chủ yếu ở giá trị rất thấp gần 0 MJ/m².

  • Temperature: Dữ liệu nhiệt độ phân bố từ -17°C tới 40°C. Trong đó, nhiệt độ phân bố chủ yếu từ 0°C tới 25°C. Giá trị 20°C có mật độ cao nhất.

  • Visibility: Biểu đồ cho thấy phần lớn các giá trị tầm nhìn tập trung ở khoảng 2000 đơn vị. Ngoài một số ít trường hợp, dữ liệu cho thấy tần suất xuất hiện của các giá trị tầm nhìn thấp hơn (dưới 2000 đơn vị) là rất ít.

  • Wind speed: Biểu đồ cho thấy tốc độ gió thường tập trung nhiều ở khoảng từ 0 đến 2m/s, với đỉnh tại khoảng 1m/s.

ggplot(data_m,aes(x = value,color = seasons)) + 
    facet_wrap(~variable, scales = "free") + 
    geom_boxplot() + theme_minimal()

Nhận xét về các boxplot:

  • Dew point temperature: Các điểm ngưng sương cao nhất ở mùa hè, thấp nhất ở mùa đông, các mùa còn lại có trung vị là trung bình của 2 mùa hè và đông

  • Humidity percent: Thấp nhất tại mùa đông và độ ẩm cao tương đương nhau ở các mùa còn lại

  • Rainfall: Lưu lượng mưa có giá trị có thể quan sát cao nhất tại mùa hè, xuân, thu còn mùa đông lượng mưa không có quan sát

  • Snowfall: Độ dày của tuyết chỉ có thể quan sát tại mùa đông

  • Solar radiation: Bức xạ mặt trời cao tại mùa hè và thấp ở mùa đông

  • Temperature: Nhiệt độ cao tại mùa hè và thấp tại mùa đông, 2 mùa còn lại tương đương nhau

  • Visbility: Cao nhất tại mùa thu và thấp nhất tại mùa xuân, mùa đông và hè tương đương

  • Wind speed: Tốc độ gió gần như không có khác biệt đáng kể giữa các mùa.

data_by_day <- data %>% mutate(date = as.Date(date, "%d/%m/%Y")) %>% group_by(date) %>% summarise(rented_bike_count = sum(rented_bike_count), seasons = first(seasons))
ggplot(data = data_by_day, aes(x=date, y=rented_bike_count, colour = seasons)) + geom_line(alpha=.3) + geom_smooth(span = 0.2, se=F) + theme_minimal()

Nhận xét:

  • Vào mùa đông, số lượng thuê xe đạp là thấp nhất.

  • Số lượng thuê xe đạp cao nhất ở thời điểm cuối mùa xuân và đầu mùa hè.

  • Số lượng thuê xe đạp có sự tăng, giảm không ổn định. Điều này có thể do tác động của những yếu tố như ngày nghỉ, lượng mưa, độ ẩm,….

ggplot(data_m,aes(y = rented_bike_count,x = value, color = time)) + 
    facet_wrap(~variable, scales = "free") + 
     geom_point( alpha = 0.3, shape = 16)+
  stat_smooth(se=F) 

Nhận xét: Nhìn chung, mối quan hệ giữa các biến độc lập và biến phụ thuộc tương đối ổn định qua các thời điểm trong ngày (trừ biến Solar radition).

ggplot(data_m,aes(y = rented_bike_count,x = value, color = seasons)) + 
    facet_wrap(~variable, scales = "free") + 
     geom_point( alpha = 0.3, shape = 16)+
  stat_smooth(se=F)

Nhận xét: Có thể thấy, mối quan hệ giữa các biến độc lập và biến phụ thuộc có sự thay đổi khá nhiều qua từng mùa.

** Vì vậy, chúng ta cần xây dựng mô hình cho từng mùa để tránh sự biến động của dữ liệu, tăng sự ổn định của mô hình hồi quy và giảm độ phức tạp trong phân tích. **

** Ngoài ra, dựa vào các biểu đồ, có thể thấy mối quan hệ giữa các biến độc lập và biến phụ tuân theo hàm mũ, do đó, chúng ta nên áp dụng phương pháp hồi quy poisson. **

ggplot(data,aes(y = humidity_percent,x = hour, color = rented_bike_count)) +
  geom_jitter(alpha = 0.5, shape = 16) 

Nhận xét:

  • Từ 0h đến 9h và từ 19h đến 23h: Có thể thấy trong hai khung giờ này, khi độ ẩm nằm trong khoảng từ 50% đến 80%, số lượng khách hàng thuê xe đạp nhiều hơn. Điều này có thể là do trong hai khung giờ đó, nhiệt độ thường khá mát mẻ hoặc lạnh, và độ ẩm từ 50% đến 80% là điều kiện lý tưởng để đi xe đạp.

  • Từ 9h đến 19h: Ở khung giờ này, thời tiết đang nóng dần lên nên khách hàng có nhu cầu thuê xe đạp nhiều hơn khi độ ẩm nằm trong khoảng từ 25% tới 60%.

ggplot(data,aes(y = temperature_c,x = hour, color = rented_bike_count)) +
  geom_jitter(alpha = 0.5, shape = 16) 

Nhận xét:

  • Số lượng xe đạp được thuê có xu hướng tăng vào các giờ sáng sớm (khoảng 7-9 giờ) và buổi chiều tối (khoảng 17-19 giờ).

  • Khi nhiệt độ tăng, đặc biệt là từ khoảng 10°C đến 30°C, số lượng xe đạp được thuê cũng tăng lên.

  • Ở các nhiệt độ cực đoan (dưới 0°C và trên 30°C), số lượng xe đạp được thuê giảm đi, điều này có thể do điều kiện thời tiết không thuận lợi cho việc đi xe đạp.

ggplot(data,aes(x = rented_bike_count,y = time)) + 
    stat_halfeye(  position = "dodge"  )+
  facet_wrap(~seasons)

Nhận xét:

  • Mùa thu: Không có sự khác biệt rõ rệt giữa các thời điểm trong ngày. Số lượng xe đạp được thuê thấp và phân bố tương đối đồng đều.

  • Mùa xuân: Buổi sáng và buổi tối có số lượng xe đạp được thuê thấp, với phân bố khá nhỏ. Buổi chiều có số lượng xe đạp được thuê cao nhất, điều này có thể do người dân tham gia nhiều hoạt động ngoài trời hơn sau giờ làm việc và khi thời tiết ấm dần lên trong ngày. Buổi trưa có số lượng xe đạp được thuê cao hơn so với sáng và tối, nhưng vẫn thấp hơn buổi chiều.

  • Mùa hè: Buổi trưa có số lượng xe đạp được thuê cao nhất. Buổi tối và buổi chiều có số lượng xe đạp được thuê thấp hơn, có thể do nhiệt độ buổi chiều và buổi tối không thuận lợi hoặc người dân ít có nhu cầu thuê xe vào các thời điểm này.

  • Mùa đông: Số lượng xe đạp được thuê rất ít vào tất cả các thời điểm trong ngày, có thể do điều kiện thời tiết khắc nghiệt, lạnh giá làm giảm nhu cầu sử dụng xe đạp.

Nhận xét: Số lượng thuê xe đạp thay đổi theo từng khung thời gian.

Nhận xét:

  • Vào mùa Thu và Xuân, phân phối số lượng xe đạp thuê có mô hình tương tự với mật độ cao hơn ở các giá trị thấp và đuôi dài mở rộng đến các giá trị cao hơn.

  • Mùa Hạ có phân phối tương đối hẹp hơn, cho thấy sự biến đổi ít hơn trong việc thuê xe đạp, với số lượng thuê vào các ngày nghỉ lễ cao hơn một chút so với các ngày không phải nghỉ lễ.

  • Mùa Đông cho thấy số lượng xe đạp thuê thấp nhất với sự giảm rõ rệt so với các mùa khác. Phân phối cho các ngày nghỉ lễ và không nghỉ lễ khá giống nhau vào mùa Đông.

Preprocessing

data.clean1 <- data |> group_by(seasons,hour) |>
  mutate(l = quantile(rented_bike_count,0.25) - 1.5*IQR(rented_bike_count), 
         h = quantile(rented_bike_count,0.75) + 1.5*IQR(rented_bike_count)) |> 
  filter(rented_bike_count >= l & rented_bike_count <= h) |> select(-c(h,l)) |> ungroup()
dis_mahalanobis <- function(X){
  return(mahalanobis(X, colMeans(X),diag(1e-9, ncol(X))+ cov(X)))
}

index <- data.clean1 |> select(-rented_bike_count)|>  mutate(id= row_number()) |> 
  group_by(seasons,hour) |> 
  mutate(mahalanobis_dist = dis_mahalanobis(across(where(is.numeric)))  ) |> 
  filter( mahalanobis_dist <= qchisq(0.95, ncol(across(where(is.numeric))))) |> 
   pull(id)
data.clean <- data.clean1[index,] |> select(-c(date))
nrow(data.clean)
## [1] 7872

Testing

Nhận định ban đầu

Khi ta xét đến các bài toán về việc quản lý các dịch vụ công, ta cần quan tâm tới những ngày mà số lượng người sử dụng dịch vụ tăng mạnh. Thông thường ta có thể thấy vào những ngày nghỉ, ngày lễ số lượng người ra đường vui chơi sẽ nhiều hơn so với ngày thường. Từ giả định đó ta sẽ kiểm tra xem việc ngày đó có phải ngày nghỉ hay không liệu có ảnh hưởng đến số lượng xe đạp được thuê hay không.

AB test cho holiday

Ta xem xét câu hỏi: “Liệu số lượng xe đạp được thuê vào các ngày nghỉ, ngày lễ có tăng đột biến hay không để từ đó thành phố có thể đề ra các kế hoạch ứng biến vào các ngày đặc biệt này”.

Trước hết ta có thể thấy rõ rằng phần trăm ngày nghỉ so với các ngày trong năm là rất thấp (ví dụ trong 1 tháng 4 tuần chỉ có 4 ngày Chủ Nhật so với 24 ngày bình thường), tuy nhiên không thể bỏ qua các trường hợp người ta sẽ gia tăng đột ngột việc thuê xe để giải trí, vui chơi vào các ngày nghỉ và ngày lễ.

Do đó ta cần kiểm định cho hai loại ngày này.

data <- data.clean
data$holiday <- as.factor(data$holiday)
data |> group_by(holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))
## # A tibble: 2 × 4
##   holiday     n  mean    sd
##   <fct>   <int> <dbl> <dbl>
## 1 0        7485  763.  643.
## 2 1         387  533.  562.

Kết quả cho thấy có 7485 lượt thuê xe vào ngày thường và 387 lượt thuê xe vào ngày nghỉ. Trung bình lượng thuê xe với ngày thường là cao hơn ngày nghỉ. Độ lệch chuẩn của ngày thường cũng lớn hơn, cho thấy độ biến động lượng thuê xe giữa các ngày thường là lớn hơn so với ngày nghỉ

Biểu đồ violin dưới đây giúp ta khẳng định các nhận định trên, đồng thời cung cấp thêm thông tin về phân phối lượng thuê xe đạp, ở đây, dữ liệu của cả hai ngày đều cho thấy phân phối bất đối xứng và lệch phải của lượng thuê xe.

ggplot(data, aes(x = holiday, y = rented_bike_count, fill = holiday)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Holiday", y = "Rented bike count") +
  theme_bw() +
  theme(legend.position = "none")

Thông qua bảng tổng hợp và biểu đổ violin, một giả định có thể được đưa ra: “số lượng thuê xe của ngày thường là nhiều hơn ngày nghỉ”. Do đó, ta cần kiểm chứng giả thuyết và đối thuyết sau:

\(H_0: \mu_1 = \mu_0\)

\(H_1: \mu_1 < \mu_0\)

Nếu \(H_0\) đúng tức sự khác biệt giữa số lượng xe đạp được thuê giữa ngày nghỉ và ngày thường chỉ là kết quả của sự ngẫu nhiên (không có ý nghĩa thống kê). Để kiểm định giả thuyết này ta sử dụng kiểm định hoán vị, và p-value sẽ được tính cho kiểm định bên trái.

perm_fun <- function(x, n1, n0, R) {
n <- n1 + n0
mean_diff <- numeric(R)
for (i in 1:R){
idx_1 <- sample(x = 1:n, size = n1)
idx_0 <- setdiff(x = 1:n, y = idx_1)
mean_diff[i] <- mean(x[idx_1]) - mean(x[idx_0])
}
return(mean_diff)
}

set.seed(42)
diff_mean_perm <- perm_fun(data$rented_bike_count, n1 = 387, n0 = 7485, R = 10000)


ggplot(data = tibble(perm_diffs = diff_mean_perm), aes(x = perm_diffs)) +
  geom_histogram(bins = 10, fill = "gray80", color = "black") +
  labs(x = "Rented bike count differences", y = "Frequency") +
  theme_bw()

Giá trị p-value

mean_0 <- mean(data$rented_bike_count[data$holiday == '0'])
mean_1 <- mean(data$rented_bike_count[data$holiday == '1'])
mean(diff_mean_perm < (mean_1 - mean_0))
## [1] 0

Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê trong ngày nghỉ và ngày thường không phải do sự ngẫu nhiên (có ý nghĩa thống kê).

Từ kết luận trên ta có thể thấy rằng, số lượng xe đạp được thuê vào ngày nghỉ sẽ không ảnh hưởng nhiều (không tăng đột biến) và vẫn ít hơn ngày thường. Tuy nhiên đây không phải yếu tố ảnh hưởng duy nhất đến việc thuê xe đạp của người dân

ANOVA test cho seasons

Tiếp đến, cùng câu hỏi như trên nhưng lần này ta lại xét đến các mùa trong năm: “số lượng thuê xe đạp có bị ảnh hưởng bởi các mùa hay không?”

Việc các mùa trong năm có thể ảnh hưởng đến số lượng thuê xe đạp rất nhiều. Ta có thể phỏng đoán rằng vào mùa đông thời tiết giá lạnh, người ta sẽ lựa chọn các loại phương tiện khác để giữ ấm thay vì chọn đạp xe dưới thời tiết khắc nghiệt kể trên. Hoặc vào mùa hè một lượng lớn du khách nước ngoài đến du lịch và trải nghiệm đạp xe ở Seoul cũng góp phần làm thay đổi số lượng thuê xe so với các mùa khác.

Do đó ta sẽ kiểm định cho cả bốn mùa trong năm.

data |> group_by(seasons) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))
## # A tibble: 4 × 4
##   seasons     n  mean    sd
##   <fct>   <int> <dbl> <dbl>
## 1 Autumn   1756  971.  604.
## 2 Spring   2057  755.  612.
## 3 Summer   2023 1085.  673.
## 4 Winter   2036  228.  148.

Đúng như chúng ta dự đoán trung bình số lượng xe đạp được thuê vào mùa đông là thấp nhất và vào mùa hè là cao nhất. Ta có thể thấy trung bình số lượng xe đạp được thuê của 4 mùa là khác nhau (về mặt giá trị).

Trong khi đó, độ lệch chuẩn là cho thấy độ biến động trong số lượng xe đạp được thuê là khác biệt giữa các nhóm.

Biểu đồ violin làm rõ hơn cho nhận xét ở trên.

ggplot(data, aes(x = seasons, y = rented_bike_count, fill = seasons)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Seasons", y = "Rented bike count") +
  theme_bw() +
  theme(legend.position = "none")

Ta đặt giả thuyết như sau:

\(H_0: \mu_{\text{autumn}} = \mu_{\text{spring}} = \mu_{\text{summer}} = \mu_{\text{winter}}\)

\(H_1:\) Ít nhất có một trung bình là khác với những cái còn lại.

Để kiểm định cho giả thuyết trên ta thực hiện Permutation ANOVA. Permutation ANOVA là một phương pháp phân tích ANOVA khác không phụ thuộc vào một số giả định như là: dữ liệu trong mỗi nhóm phải tuân theo phân phối chuẩn, giả định về sự đồng nhất phương sai,…

set.seed(42)
out_aov_1 <- aovp(formula = rented_bike_count ~ seasons, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(out_aov_1)
## Component 1 :
##               Df   R Sum Sq R Mean Sq Iter  Pr(Prob)    
## seasons        3  867441627 289147209 5000 < 2.2e-16 ***
## Residuals   7868 2369546013    301162                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Giá trị p-value được cho cung cấp bởi cột Pr(Prob) là 2.2e-16. Số lần lấy mẫu lặp lại là 5000, được cung cấp bởi cột Iter.

Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng mùa là có ý nghĩa thống kê (có ít nhất một trung bình là khác với các cái còn lại).

Từ kết luận trên, ta có thể thấy các mùa trong năm có ảnh hưởng đến số lượng thuê xe đạp. Thành phố có thể đưa ra các kế hoạch như là tăng cường xe đạp vào mùa hè để đáp ứng đủ nhu cầu và giảm bớt số lượng xe đạp vào mùa đông (có thể đưa về kho để tránh tình trạng hư hại do thời tiết).

ANOVA test kết hợp

Từ các dữ kiện ở trên ta có thể thực hiện ANOVA cho cả mùa và ngày nghỉ. Ta sẽ xem xét thêm rằng các ngày nghỉ trong một mùa có ảnh hưởng tới số lượng thuê xe đạp hay không? Ví dụ vào mùa hè các học sinh được nghỉ học tuy nhiên bố mẹ chúng vẫn đi làm, chỉ cho đến cuối tuần mới có thời gian để đưa con của họ đi chơi, số lượng thuê xe đạp có thể sẽ tăng vào những ngày như vậy.

data |> group_by(seasons, holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))
## # A tibble: 8 × 5
## # Groups:   seasons [4]
##   seasons holiday     n  mean    sd
##   <fct>   <fct>   <int> <dbl> <dbl>
## 1 Autumn  0        1663  975.  605.
## 2 Autumn  1          93  914.  581.
## 3 Spring  0        1995  758.  612.
## 4 Spring  1          62  680.  605.
## 5 Summer  0        1975 1087.  675.
## 6 Summer  1          48 1022.  564.
## 7 Winter  0        1852  235.  150.
## 8 Winter  1         184  163.  106.

Có vẻ như phỏng đoán của chúng ta lại chính xác khi mà số ngày nghỉ, ngày lễ trong mùa hè rất thấp so với ngày thường tuy nhiên trung bình lượng thuê xe đạp lại sấp sỉ với ngày thường. Điều này cũng đúng với trường hợp mùa thu (có thể do mùa này không khí dễ chịu, cảnh sắc tươi đẹp nên người ta thường dành cuối tuần để đạp xe thư giãn)

ggplot(data, aes(x = seasons, y = rented_bike_count, fill = holiday)) +
  geom_violin() +
  labs(x = "Seasons", y = "Rented bike count", fill = "Holiday") +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_wrap(~seasons, scales = "free") +
  geom_boxplot(width = 0.15, position = position_dodge(width = 0.89)) 

Để kiểm tra xem việc này là do ngẫu nhiên hay không, ta sẽ thực hiện Permutation ANOVA.

set.seed(42)
out_aov_2 <- aovp(formula = rented_bike_count ~ seasons*holiday, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(out_aov_2)
## Component 1 :
##                   Df   R Sum Sq R Mean Sq Iter Pr(Prob)    
## seasons            3  190825671  63608557 5000   <2e-16 ***
## holiday            1    1371580   1371580 5000   0.0126 *  
## seasons:holiday    3      12797      4266   51   1.0000    
## Residuals       7864 2367786918    301092                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta có thể thấy rằng, như kết luận ở phần trước, khi đứng độc lập thì các mùa trong năm và loại ngày sẽ ảnh hưởng đến số lượng thuê xe đạp (với mức ý nghĩa 5%).

Ở kiểm định này ta cũng có thể thấy thêm được sự tương tác giữa ngày nghỉ và các mùa trong năm. Với mức ý nghĩa \(\alpha = 0.05\), kết quả (p-value = 1) cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào các loại ngày theo từng mùa là không có ý nghĩa thống kê.

Do đó các ngày nghỉ, ngày lễ vẫn không phải là tác nhân chính ảnh hưởng đến việc thuê xe của người dân Hàn Quốc.

Mở rộng

Việc ngày nghỉ, ngày lễ chỉ chiếm một phần nhỏ trong năm do đó khó mà gây ảnh hưởng đến số lượng xe đạp của người dân. Ta tiếp tục xét đến yếu tố thời gian trong ngày, ta sẽ kiểm định xem thời gian trong ngày có ảnh hưởng đến số lượng thuê xe hay không

data$holiday <- as.factor(data$time)
data |> group_by(time) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))
## # A tibble: 4 × 4
##   time      n  mean    sd
##   <chr> <int> <dbl> <dbl>
## 1 chieu  1626  927.  565.
## 2 sang   2311  510.  503.
## 3 toi    3282  845.  746.
## 4 trua    653  707.  375.

Ta có thể thấy trung bình số lượng xe đạp được thuê vào buổi chiều (từ 13h đến 18h) là cao nhất

ggplot(data, aes(x = time, y = rented_bike_count, fill = time)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Seasons", y = "Rented bike count") +
  theme_bw() +
  theme(legend.position = "none")

set.seed(42)
out_aov_3 <- aovp(formula = rented_bike_count ~ time, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(out_aov_3)
## Component 1 :
##               Df   R Sum Sq R Mean Sq Iter  Pr(Prob)    
## time1          3  214695712  71565237 5000 < 2.2e-16 ***
## Residuals   7868 3022291929    384125                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng khung giờ là có ý nghĩa thống kê (có ít nhất một trung bình là khác với các cái còn lại).

Từ kết luận trên, ta có thể thấy các khung giờ trong ngày có ảnh hưởng đến số lượng thuê xe đạp. Thành phố có thể đưa ra các kế hoạch như là tăng cường xe đạp vào từng khung giờ cụ thể để đáp ứng đủ nhu cầu của người dân.

Tương tự như trên ta cũng kiểm tra xem các khung giờ theo từng mùa trong năm

data |> group_by(seasons, holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))
## # A tibble: 16 × 5
## # Groups:   seasons [4]
##    seasons holiday     n  mean    sd
##    <fct>   <fct>   <int> <dbl> <dbl>
##  1 Autumn  chieu     356 1276. 421. 
##  2 Autumn  sang      512  661. 548. 
##  3 Autumn  toi       750 1038. 665. 
##  4 Autumn  trua      138  974. 212. 
##  5 Spring  chieu     429 1044. 541. 
##  6 Spring  sang      594  491. 478. 
##  7 Spring  toi       862  789. 689. 
##  8 Spring  trua      172  777. 357. 
##  9 Summer  chieu     428 1120. 488. 
## 10 Summer  sang      599  730. 536. 
## 11 Summer  toi       821 1374. 762. 
## 12 Summer  trua      175  865. 275. 
## 13 Winter  chieu     413  306. 105. 
## 14 Winter  sang      606  184. 171. 
## 15 Winter  toi       849  218. 144. 
## 16 Winter  trua      168  253.  78.7

Tất cả các mùa đều có trung bình số lượng xe đạp được thuê vào buổi chiều là cao nhất.

ggplot(data, aes(x = seasons, y = rented_bike_count, fill = time)) +
  geom_violin() +
  labs(x = "Seasons", y = "Rented bike count", fill = "Time") +
  theme_bw() +
  theme(legend.position = "bottom") +
  facet_wrap(~seasons, scales = "free") +
  geom_boxplot(width = 0.15, position = position_dodge(width = 0.89)) 

set.seed(42)
out_aov_4 <- aovp(formula = rented_bike_count ~ seasons*time, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(out_aov_4)
## Component 1 :
##                 Df   R Sum Sq R Mean Sq Iter  Pr(Prob)    
## seasons          3  548011049 182670350 5000 < 2.2e-16 ***
## time1            3  217347833  72449278 5000 < 2.2e-16 ***
## seasons:time1    9  109074873  12119430 5000 < 2.2e-16 ***
## Residuals     7856 2048677419    260779                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng khung giờ theo từng mùa là có ý nghĩa thống kê, cho thấy rằng tác động của mùa lên số lượng xe đạp thuê thay đổi tùy thuộc vào từng khung giờ.

Regression

Data Model

Xây dựng mô hình hồi quy tuyến tính với biến phụ thuộc là rented_bike_count, và biến giải thích là các biến còn lại. 2 hướng nhắm tới là Linear Models và Generalized Linear Models. Chia tập dữ liệu thành 2 phần để đánh giá.

data.model <- data.clean |> select(-time)

slit_test_train <- function(df,pro=0.8,seed =30){
  set.seed(seed)
  N <-  nrow(df)
  ind_train <- sample(1:N,size=floor(pro*N))
  ind_test <- setdiff(1:N,ind_train)
  data_train <- tibble( df[ind_train,])
  data_test <- tibble( df[ind_test,])
  return(list(data_train,data_test))
}
dummies <- dummyVars("~ .", data = data.model |> select(-rented_bike_count))
data.numeric <- data.frame(predict(dummies, newdata = data.model |> select(-rented_bike_count)))
data.numeric$rented_bike_count <- data.model$rented_bike_count
rs <- slit_test_train(data.numeric,0.8)
trainset_numeric <- as.data.frame(rs[1])
testset_numeric <- as.data.frame(rs[2])
rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])

Linear Regression

md.ln <- lm(rented_bike_count ~ .,
                           data = trainset)
summary(md.ln)
## 
## Call:
## lm(formula = rented_bike_count ~ ., data = trainset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1352.97  -215.39   -11.04   200.77  1725.97 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              1.306e+03  1.230e+02  10.618  < 2e-16 ***
## hour1                   -1.231e+02  3.211e+01  -3.835 0.000127 ***
## hour2                   -2.399e+02  3.199e+01  -7.499 7.31e-14 ***
## hour3                   -3.060e+02  3.191e+01  -9.589  < 2e-16 ***
## hour4                   -3.816e+02  3.182e+01 -11.993  < 2e-16 ***
## hour5                   -3.580e+02  3.229e+01 -11.084  < 2e-16 ***
## hour6                   -1.970e+02  3.217e+01  -6.126 9.58e-10 ***
## hour7                    1.423e+02  3.211e+01   4.431 9.54e-06 ***
## hour8                    4.603e+02  3.228e+01  14.260  < 2e-16 ***
## hour9                    3.937e+01  3.319e+01   1.186 0.235628    
## hour10                  -1.783e+02  3.488e+01  -5.111 3.30e-07 ***
## hour11                  -1.888e+02  3.606e+01  -5.235 1.70e-07 ***
## hour12                  -1.321e+02  3.691e+01  -3.580 0.000346 ***
## hour13                  -1.439e+02  3.776e+01  -3.813 0.000139 ***
## hour14                  -1.487e+02  3.637e+01  -4.088 4.41e-05 ***
## hour15                  -5.469e+01  3.586e+01  -1.525 0.127300    
## hour16                   1.149e+02  3.464e+01   3.318 0.000911 ***
## hour17                   3.633e+02  3.347e+01  10.855  < 2e-16 ***
## hour18                   8.332e+02  3.277e+01  25.427  < 2e-16 ***
## hour19                   5.633e+02  3.269e+01  17.233  < 2e-16 ***
## hour20                   4.896e+02  3.239e+01  15.115  < 2e-16 ***
## hour21                   4.643e+02  3.225e+01  14.398  < 2e-16 ***
## hour22                   3.768e+02  3.208e+01  11.746  < 2e-16 ***
## hour23                   1.262e+02  3.186e+01   3.959 7.60e-05 ***
## temperature_c            3.806e+00  4.736e+00   0.804 0.421647    
## humidity_percent        -1.010e+01  1.347e+00  -7.499 7.32e-14 ***
## wind_speed_m_s           5.202e+00  5.424e+00   0.959 0.337584    
## visibility_10m          -3.460e-03  9.899e-03  -0.350 0.726679    
## dew_point_temperature_c  2.007e+01  4.986e+00   4.025 5.76e-05 ***
## solar_radiation_mj_m2    6.831e+01  1.147e+01   5.955 2.74e-09 ***
## rainfall_mm             -2.526e+02  1.956e+01 -12.917  < 2e-16 ***
## snowfall_cm              6.868e+01  1.945e+01   3.531 0.000418 ***
## seasonsSpring           -1.824e+02  1.383e+01 -13.188  < 2e-16 ***
## seasonsSummer           -1.750e+02  1.724e+01 -10.154  < 2e-16 ***
## seasonsWinter           -4.089e+02  1.969e+01 -20.770  < 2e-16 ***
## holiday1                -1.356e+02  2.120e+01  -6.396 1.71e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 358.5 on 6261 degrees of freedom
## Multiple R-squared:  0.6835, Adjusted R-squared:  0.6818 
## F-statistic: 386.4 on 35 and 6261 DF,  p-value: < 2.2e-16

Mô hình giải thích được hơn 68% lượng thay đổi trong các biến giải thích

anova(md.ln, test="Chisq")
## Analysis of Variance Table
## 
## Response: rented_bike_count
##                           Df    Sum Sq   Mean Sq   F value    Pr(>F)    
## hour                      23 852915014  37083261  288.5345 < 2.2e-16 ***
## temperature_c              1 699768135 699768135 5444.7008 < 2.2e-16 ***
## humidity_percent           1  62859972  62859972  489.0959 < 2.2e-16 ***
## wind_speed_m_s             1     83726     83726    0.6514 0.4196257    
## visibility_10m             1   1030162   1030162    8.0154 0.0046529 ** 
## dew_point_temperature_c    1   1934492   1934492   15.0517 0.0001057 ***
## solar_radiation_mj_m2      1   6674927   6674927   51.9357 6.409e-13 ***
## rainfall_mm                1  22859938  22859938  177.8668 < 2.2e-16 ***
## snowfall_cm                1      3532      3532    0.0275 0.8683299    
## seasons                    3  84629981  28209994  219.4941 < 2.2e-16 ***
## holiday                    1   5257906   5257906   40.9103 1.709e-10 ***
## Residuals               6261 804681192    128523                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Các biến hour, temperature_c, humidity_percent, dew_point_temperature_c, solar_radiation_mj_m2, rainfall_mm, seasons, và holiday có ý nghĩa thống kê cao và ảnh hưởng lớn đến số lượng xe đạp thuê.
  • Các biến wind_speed_m_s và snowfall_cm không có ý nghĩa thống kê trong mô hình này.
y_pre_train <- predict(md.ln,trainset,type = "response")
y_pre_test <- predict(md.ln,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("RMSE của trainset:",rmse_train))
## [1] "RMSE của trainset: 357.474509243905"
print(paste("RMSE của testset:",rmse_test))
## [1] "RMSE của testset: 369.427222345817"
md.ln.cv <- errorest(rented_bike_count ~ .,
                         data = data.model,model = lm,estimator="cv",
                     est.para=control.errorest(k=5, predictions = TRUE))
md.ln.cv
## 
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model, 
##     model = lm, estimator = "cv", est.para = control.errorest(k = 5, 
##         predictions = TRUE))
## 
##   5-fold cross-validation estimator of root mean squared error
## 
## Root mean squared error:  361.5928

Dựa vào Rmse trên test, train, cũng như kết quả k-fold cross-validation, thì kết quả mô hình không bị overfitting nhưng rmse vẫn còn cao.

set.seed(12)
rmse_l <- rep(1,ncol(data.numeric))
for( i in 1:ncol(data.numeric)){
  data.pca <- as.data.frame(PCA(data.numeric |> select(-rented_bike_count), ncp =i, graph = F)$ind$coord) 
  data.pca$rented_bike_count <- data.numeric$rented_bike_count
  md.ln.cv <- errorest(rented_bike_count ~ .,
                           data = data.pca,model = lm,estimator="cv",
                       est.para=control.errorest(k=5, predictions = TRUE))
  rmse_l[i] <- md.ln.cv$error
}
which.min(rmse_l)
## [1] 38
rmse_l
##  [1] 572.9073 522.7233 522.5773 509.0930 498.7240 461.3164 447.7140 447.0399
##  [9] 420.6869 418.3386 416.4121 416.5305 415.0307 413.1027 412.1735 410.1889
## [17] 407.2417 402.2952 401.9710 402.2343 400.7602 401.0530 400.9924 400.4795
## [25] 400.7529 397.9674 397.5725 397.8474 376.2425 376.5101 367.1026 363.5321
## [33] 361.7508 361.5474 361.6642 361.4610 361.4258 361.4140 361.5543

Dù có PCA thì kết quả vẫn không hiệu quả hơn là mấy, vậy có thể kết luận dữ liệu không bị đa công tuyến quá nhiều.

Poisson regression

md.ps <- glm(rented_bike_count ~ .,
                           data = trainset,family = poisson)
summary(md.ps)
## 
## Call:
## glm(formula = rented_bike_count ~ ., family = poisson, data = trainset)
## 
## Coefficients:
##                           Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)              7.301e+00  1.337e-02  546.144  < 2e-16 ***
## hour1                   -2.439e-01  3.965e-03  -61.524  < 2e-16 ***
## hour2                   -5.692e-01  4.364e-03 -130.405  < 2e-16 ***
## hour3                   -9.130e-01  4.967e-03 -183.813  < 2e-16 ***
## hour4                   -1.360e+00  5.805e-03 -234.262  < 2e-16 ***
## hour5                   -1.283e+00  5.814e-03 -220.752  < 2e-16 ***
## hour6                   -5.431e-01  4.453e-03 -121.955  < 2e-16 ***
## hour7                    2.108e-01  3.615e-03   58.324  < 2e-16 ***
## hour8                    5.965e-01  3.348e-03  178.167  < 2e-16 ***
## hour9                    1.125e-01  3.728e-03   30.174  < 2e-16 ***
## hour10                  -1.572e-01  4.063e-03  -38.682  < 2e-16 ***
## hour11                  -1.179e-01  4.124e-03  -28.587  < 2e-16 ***
## hour12                  -4.731e-03  4.107e-03   -1.152  0.24933    
## hour13                  -1.286e-02  4.156e-03   -3.095  0.00197 ** 
## hour14                  -1.752e-02  4.029e-03   -4.348 1.37e-05 ***
## hour15                   8.817e-02  3.894e-03   22.645  < 2e-16 ***
## hour16                   2.524e-01  3.652e-03   69.119  < 2e-16 ***
## hour17                   4.728e-01  3.424e-03  138.077  < 2e-16 ***
## hour18                   8.230e-01  3.200e-03  257.208  < 2e-16 ***
## hour19                   6.486e-01  3.257e-03  199.117  < 2e-16 ***
## hour20                   5.929e-01  3.272e-03  181.198  < 2e-16 ***
## hour21                   5.860e-01  3.283e-03  178.504  < 2e-16 ***
## hour22                   4.884e-01  3.315e-03  147.314  < 2e-16 ***
## hour23                   1.975e-01  3.523e-03   56.046  < 2e-16 ***
## temperature_c            1.589e-03  4.956e-04    3.206  0.00134 ** 
## humidity_percent        -1.223e-02  1.496e-04  -81.762  < 2e-16 ***
## wind_speed_m_s          -1.028e-02  5.772e-04  -17.805  < 2e-16 ***
## visibility_10m          -2.278e-05  1.044e-06  -21.820  < 2e-16 ***
## dew_point_temperature_c  2.635e-02  5.240e-04   50.292  < 2e-16 ***
## solar_radiation_mj_m2    3.165e-02  1.171e-03   27.021  < 2e-16 ***
## rainfall_mm             -1.246e+00  6.741e-03 -184.888  < 2e-16 ***
## snowfall_cm             -7.930e-02  4.059e-03  -19.537  < 2e-16 ***
## seasonsSpring           -2.068e-01  1.351e-03 -153.076  < 2e-16 ***
## seasonsSummer           -2.033e-01  1.601e-03 -126.974  < 2e-16 ***
## seasonsWinter           -1.019e+00  2.433e-03 -418.852  < 2e-16 ***
## holiday1                -1.812e-01  2.545e-03  -71.194  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 3247619  on 6296  degrees of freedom
## Residual deviance:  640109  on 6261  degrees of freedom
## AIC: 690550
## 
## Number of Fisher Scoring iterations: 5

Dựa vào bảng thống kê trên ta nhận thấy tham số phân tán là 1, nhưng thực chất với dữ liệu thực tế tình trạng này là thường xuyên gặp. Cụ thể trong trường hợp này là hơn 100 (\(\phi = \frac{\sum (\text{Pearson residuals})^2}{n - p}\)), do đấy mô hình đang bị overdispersion. Khi đó giả định về giả định rằng phương sai bằng trung bình không còn đúng.

anova(md.ps, test="Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: rented_bike_count
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                     6296    3247619              
## hour                    23  1248049      6273    1999569 < 2.2e-16 ***
## temperature_c            1   963614      6272    1035955 < 2.2e-16 ***
## humidity_percent         1    46470      6271     989485 < 2.2e-16 ***
## wind_speed_m_s           1     2585      6270     986900 < 2.2e-16 ***
## visibility_10m           1      234      6269     986666 < 2.2e-16 ***
## dew_point_temperature_c  1     4743      6268     981923 < 2.2e-16 ***
## solar_radiation_mj_m2    1     1657      6267     980266 < 2.2e-16 ***
## rainfall_mm              1    72578      6266     907688 < 2.2e-16 ***
## snowfall_cm              1    12412      6265     895276 < 2.2e-16 ***
## seasons                  3   249812      6262     645464 < 2.2e-16 ***
## holiday                  1     5355      6261     640109 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Các biến giải thích đều có p-value < 2.2e-16.

  • Biến giờ (hour), nhiệt độ (temperature_c), độ ẩm (humidity_percent), mùa (seasons) và lượng mưa (rainfall_mm) là những biến có ảnh hưởng lớn nhất.

  • Các biến khác như tốc độ gió, tầm nhìn, nhiệt độ điểm sương, bức xạ mặt trời, lượng tuyết rơi, và ngày lễ cũng có ảnh hưởng nhưng mức độ ít hơn so với các biến chính.

y_pre_train <- predict(md.ps,trainset,type = "response")
y_pre_test <- predict(md.ps,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 312.076793691618"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 320.641281267276"
myGLM = function(formula, data) {
  glm(formula, data, family = poisson(link = log))
}

myPredictGLM = function(object, newdata){
  predict(object, newdata , type="response")
}

md.ps.cv <- errorest(rented_bike_count ~ .,
                         data = data.model, predict  = myPredictGLM,model = myGLM,estimator="cv",
                     est.para=control.errorest(k=5, predictions = TRUE))
md.ps.cv
## 
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model, 
##     model = myGLM, predict = myPredictGLM, estimator = "cv", 
##     est.para = control.errorest(k = 5, predictions = TRUE))
## 
##   5-fold cross-validation estimator of root mean squared error
## 
## Root mean squared error:  315.2939

Kết quả rmse tốt hơn so với mô hình hồi quy tuyến tính.

md.ps.null <- glm(rented_bike_count ~ 1, family = "poisson", data = trainset)
md.ps.R2 <- 1 - as.numeric(logLik(md.ps)) / as.numeric(logLik(md.ps.null))
md.ps.R2
## [1] 0.7906365

Với \(Pseudo-R^2\) gần 0.8 cho thấy mức giải thích được cải thiện là gần 80% với vo mô hình chỉ có bias.

Quasipoisson

md.qps <- glm(rented_bike_count ~ .,
                           data = trainset,family = quasipoisson)
summary(md.qps)
## 
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson, data = trainset)
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              7.301e+00  1.348e-01  54.148  < 2e-16 ***
## hour1                   -2.439e-01  3.999e-02  -6.100 1.13e-09 ***
## hour2                   -5.692e-01  4.402e-02 -12.929  < 2e-16 ***
## hour3                   -9.130e-01  5.010e-02 -18.224  < 2e-16 ***
## hour4                   -1.360e+00  5.856e-02 -23.226  < 2e-16 ***
## hour5                   -1.283e+00  5.864e-02 -21.887  < 2e-16 ***
## hour6                   -5.431e-01  4.492e-02 -12.091  < 2e-16 ***
## hour7                    2.108e-01  3.646e-02   5.783 7.71e-09 ***
## hour8                    5.965e-01  3.377e-02  17.664  < 2e-16 ***
## hour9                    1.125e-01  3.760e-02   2.992 0.002786 ** 
## hour10                  -1.572e-01  4.098e-02  -3.835 0.000127 ***
## hour11                  -1.179e-01  4.159e-02  -2.834 0.004608 ** 
## hour12                  -4.731e-03  4.142e-02  -0.114 0.909073    
## hour13                  -1.286e-02  4.192e-02  -0.307 0.758998    
## hour14                  -1.752e-02  4.064e-02  -0.431 0.666429    
## hour15                   8.817e-02  3.927e-02   2.245 0.024795 *  
## hour16                   2.524e-01  3.683e-02   6.853 7.93e-12 ***
## hour17                   4.728e-01  3.453e-02  13.690  < 2e-16 ***
## hour18                   8.230e-01  3.227e-02  25.501  < 2e-16 ***
## hour19                   6.486e-01  3.285e-02  19.741  < 2e-16 ***
## hour20                   5.929e-01  3.300e-02  17.965  < 2e-16 ***
## hour21                   5.860e-01  3.311e-02  17.698  < 2e-16 ***
## hour22                   4.884e-01  3.344e-02  14.605  < 2e-16 ***
## hour23                   1.975e-01  3.553e-02   5.557 2.86e-08 ***
## temperature_c            1.589e-03  4.999e-03   0.318 0.750575    
## humidity_percent        -1.223e-02  1.509e-03  -8.106 6.23e-16 ***
## wind_speed_m_s          -1.028e-02  5.822e-03  -1.765 0.077573 .  
## visibility_10m          -2.278e-05  1.053e-05  -2.163 0.030552 *  
## dew_point_temperature_c  2.635e-02  5.285e-03   4.986 6.32e-07 ***
## solar_radiation_mj_m2    3.165e-02  1.181e-02   2.679 0.007403 ** 
## rainfall_mm             -1.246e+00  6.800e-02 -18.331  < 2e-16 ***
## snowfall_cm             -7.930e-02  4.094e-02  -1.937 0.052786 .  
## seasonsSpring           -2.068e-01  1.362e-02 -15.177  < 2e-16 ***
## seasonsSummer           -2.033e-01  1.615e-02 -12.589  < 2e-16 ***
## seasonsWinter           -1.019e+00  2.454e-02 -41.527  < 2e-16 ***
## holiday1                -1.812e-01  2.567e-02  -7.059 1.86e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 101.7313)
## 
##     Null deviance: 3247619  on 6296  degrees of freedom
## Residual deviance:  640109  on 6261  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
anova(md.qps, test="Chisq")
## Analysis of Deviance Table
## 
## Model: quasipoisson, link: log
## 
## Response: rented_bike_count
## 
## Terms added sequentially (first to last)
## 
## 
##                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                     6296    3247619              
## hour                    23  1248049      6273    1999569 < 2.2e-16 ***
## temperature_c            1   963614      6272    1035955 < 2.2e-16 ***
## humidity_percent         1    46470      6271     989485 < 2.2e-16 ***
## wind_speed_m_s           1     2585      6270     986900 4.640e-07 ***
## visibility_10m           1      234      6269     986666    0.1293    
## dew_point_temperature_c  1     4743      6268     981923 8.591e-12 ***
## solar_radiation_mj_m2    1     1657      6267     980266 5.449e-05 ***
## rainfall_mm              1    72578      6266     907688 < 2.2e-16 ***
## snowfall_cm              1    12412      6265     895276 < 2.2e-16 ***
## seasons                  3   249812      6262     645464 < 2.2e-16 ***
## holiday                  1     5355      6261     640109 4.009e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Các biến hour, temperature_c, humidity_percent, wind_speed_m_s, dew_point_temperature_c, solar_radiation_mj_m2, rainfall_mm, snowfall_cm, seasons, và holiday đều có mức ý nghĩa thống kê cao (Pr(>Chi) < 0.05).
  • Biến visibility_10m không có mức ý nghĩa thống kê (Pr(>Chi) > 0.05), cho thấy nó không có ảnh hưởng đáng kể đến số lượng xe đạp thuê trong mô hình này.
  • Dữ liệu còn khá phân tán có thể chia ra thành các mùa, để các dữ liệu ít bị phân tác hơn.
y_pre_train <- predict(md.qps,trainset,type = "response")
y_pre_test <- predict(md.qps,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 312.076793691618"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 320.641281267276"

Đánh giá

autoplot(md.ps, which = 1, ncol = 1, label.size = 3,
         colour = "seasons",alpha = 0.6) + theme_bw()

Các điểm thặng dư đối xứng qua trục y=0, không xuất hiện đường cong bất thường nào. Điều này cho thấy tính tuyến tính của mô hình được bảo đảm. Tuy nhiên giả định đồng nhất phương sai là không đảm bảo, vì thặng dư co cụm về bên phải tăng đều, và phân tán mạnh hơn. Điều này một phần vì dữ liệu phân tán khá lớn, tuy nhiên các mô hình cũng không cần giả định này phải thoả mãn

autoplot(md.ps, which = 2, ncol = 1, label.size = 3,
         colour = "seasons",alpha = 0.6) + theme_bw()

Thặng dư không tuân theo phân phối chuẩn, điều này là hiển nhiên vì đây là hồi quy poisson, nhưng số lượng mẫu khá lớn nên thăng dư cũng xấp xỉ gần với phân phối chuẩn.

autoplot(md.ps, which = 3, ncol = 1, label.size = 3,
         colour = "seasons",alpha = 0.6) + theme_bw()

Như đã đề cập ở dữ liệu phân tán khá cao, giả định về đồng nhất phương sai thặng dư không đảm bảo. Ở đây ta thấy rõ ràng ngoài xu hướng càng tăng giá trị fit, càng tăng dần độ phân tán, độ phân tán còn có sự khác nhau ở các mùa.

autoplot(md.ps, which = 4:6, ncol = 1, label.size = 3,
         colour = "seasons",alpha = 0.6) + theme_bw()

  • Cook’s Distance: dùng để đo lường mức độ ảnh hưởng của từng điểm dữ liệu đối với mô hình hồi quy. Điểm dữ liệu có giá trị Cook’s Distance lớn có thể có ảnh hưởng lớn đến các tham số ước lượng của mô hình. +) Điểm 5212, 4494, 6029: Các điểm này có giá trị Cook’s Distance cao hơn hẳn so với các điểm khác, cho thấy chúng có ảnh hưởng lớn đến mô hình. Những điểm này có thể là outliers hoặc có leverage cao. +) Phần lớn các điểm có giá trị Cook’s Distance nhỏ, cho thấy chúng không ảnh hưởng nhiều đến mô hình.

  • Residuals vs Leverage: kiểm tra mối quan hệ giữa residuals (sai số dự đoán) và leverage (mức độ ảnh hưởng của một quan sát trong việc xác định ước lượng của mô hình). +) Residuals lớn (Điểm 5212, 6029, 4494): Các điểm này có residuals lớn, có thể là dấu hiệu của các outliers, cho thấy mô hình dự đoán kém đối với những điểm này. +) Phân bố residuals: Phân bố residuals không đối xứng hoàn toàn xung quanh trục y = 0. Điều này có thể chỉ ra rằng mô hình có vấn đề về phân phối sai số. +) Leverage thấp: Hầu hết các điểm có leverage thấp (gần 0), cho thấy chúng không ảnh hưởng mạnh đến mô hình.

  • Cook’s Distance vs Leverage: kiểm tra mối quan hệ giữa Cook’s Distance và leverage. Lệu những điểm có leverage cao (ảnh hưởng lớn) cũng có Cook’s Distance cao (ảnh hưởng đến mô hình) hay không. +) Điểm 5212, 6029, 4494: Các điểm này có cả Cook’s Distance và leverage cao, cho thấy chúng có ảnh hưởng lớn đến mô hình và cần được xem xét kỹ lưỡng. Nếu các điểm này là các outliers, có thể cân nhắc loại bỏ chúng hoặc kiểm tra kỹ hơn để xem chúng có hợp lý hay không. +) Phần lớn các điểm: Hầu hết các điểm còn lại có Cook’s Distance và leverage thấp, cho thấy chúng không ảnh hưởng nhiều đến mô hình.

data.md.vs <- data.frame(trainset |> select(temperature_c,humidity_percent,wind_speed_m_s,visibility_10m,rainfall_mm,snowfall_cm,seasons) , parital_residual = residuals(md.ps) )
data.md.vs.m <- gather(data.md.vs,var,val,-c(seasons,parital_residual) )
ggplot(data.md.vs.m,aes(x = val, y  = parital_residual, color = seasons))+facet_wrap(~var,scales = 'free_x') + geom_point()

  • Humidity Percent +) Các giá trị partial residuals phân bố khá đồng đều quanh trục y = 0. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Độ ẩm có vẻ không có mối quan hệ phi tuyến với biến đáp ứng. Có thể biến này được mô hình hóa tốt bởi Poisson Regression.

  • Rainfall mm +) Các giá trị partial residuals khá phân tán, đặc biệt ở mức lượng mưa thấp (dưới 1 mm). +) Phân bố partial residuals không thay đổi nhiều theo mùa. +) Mối quan hệ không rõ ràng, có thể do một số outliers. Cần kiểm tra thêm hoặc có thể cân nhắc biến đổi biến này

  • Snowfall cm +) Phân bố partial residuals: Các giá trị partial residuals chủ yếu tập trung ở mức tuyết rơi bằng 0, với một số outliers ở các giá trị tuyết rơi cao hơn. +) Phần lớn các điểm tuyết rơi lớn xuất hiện vào mùa đông. +) Mối quan hệ phi tuyến rõ ràng hơn, đặc biệt vào mùa đông. Cần kiểm tra thêm và có thể cần một số biến đổi hoặc thêm biến tương tác để cải thiện mô hình.

  • Temperature c +) Các giá trị partial residuals phân bố khá đồng đều quanh trục y = 0. +) Có sự phân bố rõ ràng theo mùa: nhiệt độ cao hơn vào mùa hè và thấp hơn vào mùa đông. +) Mối quan hệ giữa nhiệt độ và partial residuals có thể được biểu diễn bằng một mô hình tuyến tính.

  • Visibility 10m +) Các giá trị partial residuals phân bố đồng đều, với một số outliers. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Mối quan hệ không rõ ràng. Cần kiểm tra thêm để xem liệu biến này có cần biến đổi hoặc loại bỏ không.

  • Wind Speed m/s +) Các giá trị partial residuals phân bố đồng đều quanh trục y = 0. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Tốc độ gió có vẻ được mô hình hóa tốt bởi Poisson Regression. Mối quan hệ tuyến tính có thể hợp lý.

Khoảng tin cậy bootstraping:

fun.boot<- function(data, ind, formula, ...){
  data_new <- data[ind,]
  out_md <- glm(formula = formula, data = data_new,family = poisson,...) 
  return(out_md$coefficients)
}
md.ps.boot <- boot(data = data.model, statistic = fun.boot, R = 500,
  formula = rented_bike_count~.)
summary(md.ps.boot)
## 
## Number of bootstrap replications R = 500 
##       original    bootBias     bootSE     bootMed
## 1   7.2013e+00 -3.7976e-03 1.3985e-01  7.1927e+00
## 2  -2.4136e-01  2.8829e-03 2.5424e-02 -2.3931e-01
## 3  -5.7418e-01  1.4716e-03 2.7498e-02 -5.7214e-01
## 4  -9.2478e-01  2.4856e-04 3.0532e-02 -9.2251e-01
## 5  -1.3666e+00  1.9385e-03 2.7247e-02 -1.3637e+00
## 6  -1.2869e+00  3.2671e-04 2.8212e-02 -1.2847e+00
## 7  -5.3392e-01  2.0532e-03 3.2465e-02 -5.3123e-01
## 8   1.9199e-01  2.9494e-03 3.4654e-02  1.9531e-01
## 9   6.3545e-01  2.3324e-03 3.5198e-02  6.3773e-01
## 10  1.1536e-01  1.5840e-03 2.7225e-02  1.1463e-01
## 11 -1.7923e-01  2.5840e-03 2.9299e-02 -1.7710e-01
## 12 -1.1787e-01  2.3996e-03 3.1779e-02 -1.1614e-01
## 13 -1.2281e-02  1.8045e-03 3.4071e-02 -1.0055e-02
## 14 -2.6172e-02  1.3404e-03 3.5162e-02 -2.4951e-02
## 15 -2.1157e-02  2.7690e-03 3.3273e-02 -1.8601e-02
## 16  7.7703e-02  2.9788e-03 3.4197e-02  8.1346e-02
## 17  2.3722e-01  6.7044e-04 3.2202e-02  2.3882e-01
## 18  4.8965e-01  3.8637e-04 2.6837e-02  4.8986e-01
## 19  8.2282e-01  2.6467e-03 2.5710e-02  8.2455e-01
## 20  6.5895e-01  1.5563e-03 2.7639e-02  6.6107e-01
## 21  5.9000e-01  2.1887e-03 2.5735e-02  5.9137e-01
## 22  5.8623e-01  2.3546e-03 2.5011e-02  5.8968e-01
## 23  4.8602e-01  2.3758e-03 2.6354e-02  4.8956e-01
## 24  1.9370e-01  1.7097e-03 2.6574e-02  1.9476e-01
## 25  4.1078e-03  9.1163e-05 5.0174e-03  4.2269e-03
## 26 -1.1191e-02  3.0218e-05 1.6019e-03 -1.1167e-02
## 27 -8.9989e-03 -1.2178e-04 5.5612e-03 -9.0632e-03
## 28 -1.3289e-05 -4.8165e-08 1.0229e-05 -1.3441e-05
## 29  2.3357e-02 -8.0026e-05 5.2639e-03  2.3203e-02
## 30  3.3739e-02 -3.1836e-04 1.1250e-02  3.3410e-02
## 31 -1.3279e+00 -7.2986e-03 1.2356e-01 -1.3333e+00
## 32 -8.1502e-02 -4.2052e-04 2.0383e-02 -8.0285e-02
## 33 -1.9582e-01 -1.1999e-03 1.3003e-02 -1.9600e-01
## 34 -1.9601e-01 -1.7357e-04 1.5688e-02 -1.9461e-01
## 35 -1.0123e+00  4.2797e-04 2.2067e-02 -1.0117e+00
## 36 -1.7515e-01 -1.9252e-04 3.0740e-02 -1.7603e-01
boot.ci(md.ps.boot, index = 1, type = "perc", conf = 0.95)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 500 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = md.ps.boot, conf = 0.95, type = "perc", index = 1)
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 6.939,  7.495 )  
## Calculations and Intervals on Original Scale
set.seed(14)
sample <- data.model |> select(-rented_bike_count) |> sample_n(1) |> mutate(across(where(is.numeric), ~ . + .*rnorm(1, mean = 0, sd = 0.01)))
print(t(sample))
##                         [,1]       
## hour                    "21"       
## temperature_c           "-8.831034"
## humidity_percent        "43.0568"  
## wind_speed_m_s          "1.400403" 
## visibility_10m          "1995.394" 
## dew_point_temperature_c "-19.04248"
## solar_radiation_mj_m2   "0"        
## rainfall_mm             "0"        
## snowfall_cm             "0"        
## seasons                 "Winter"   
## holiday                 "0"
sample <- model.matrix(~ ., data = sample)
y.pre.boot <- apply(md.ps.boot$t, 1, function(x){
  exp(sample %*% x)
  }) 
quantile(y.pre.boot, probs = c(0.025, 0.975))
##     2.5%    97.5% 
## 307.7393 335.3752

Cải thiện model

Mô hình cho từng mùa

Winter

data.model.wt <- data.model[ data.model$ seasons =="Winter",- which(colnames(data.model) %in%  c('seasons'))]
rs <- slit_test_train(data.model.wt,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])


md.qps.wt <- glm(rented_bike_count ~  .,family = quasipoisson('log'),data=trainset )
summary(md.qps.wt)
## 
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson("log"), 
##     data = trainset)
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              6.102e+00  2.293e-01  26.610  < 2e-16 ***
## hour1                   -7.128e-03  6.951e-02  -0.103 0.918331    
## hour2                   -3.884e-01  7.420e-02  -5.234 1.87e-07 ***
## hour3                   -7.416e-01  8.407e-02  -8.822  < 2e-16 ***
## hour4                   -1.167e+00  9.796e-02 -11.909  < 2e-16 ***
## hour5                   -1.132e+00  9.463e-02 -11.958  < 2e-16 ***
## hour6                   -4.981e-01  7.681e-02  -6.485 1.18e-10 ***
## hour7                    3.253e-01  6.351e-02   5.121 3.41e-07 ***
## hour8                    9.490e-01  5.984e-02  15.858  < 2e-16 ***
## hour9                    5.194e-01  6.270e-02   8.284 2.50e-16 ***
## hour10                   1.298e-01  6.823e-02   1.902 0.057336 .  
## hour11                   2.351e-01  6.917e-02   3.398 0.000694 ***
## hour12                   3.132e-01  7.061e-02   4.435 9.82e-06 ***
## hour13                   3.348e-01  7.299e-02   4.587 4.85e-06 ***
## hour14                   3.292e-01  6.978e-02   4.718 2.59e-06 ***
## hour15                   3.527e-01  6.733e-02   5.238 1.83e-07 ***
## hour16                   3.965e-01  6.548e-02   6.055 1.74e-09 ***
## hour17                   5.296e-01  6.146e-02   8.618  < 2e-16 ***
## hour18                   8.481e-01  5.839e-02  14.525  < 2e-16 ***
## hour19                   4.945e-01  6.119e-02   8.081 1.26e-15 ***
## hour20                   2.730e-01  6.418e-02   4.253 2.23e-05 ***
## hour21                   3.054e-01  6.262e-02   4.878 1.18e-06 ***
## hour22                   2.630e-01  6.477e-02   4.060 5.13e-05 ***
## hour23                   3.175e-02  6.679e-02   0.475 0.634626    
## temperature_c            1.591e-02  8.364e-03   1.903 0.057282 .  
## humidity_percent        -1.102e-02  2.589e-03  -4.256 2.20e-05 ***
## wind_speed_m_s          -3.345e-02  8.443e-03  -3.962 7.77e-05 ***
## visibility_10m           1.625e-04  2.537e-05   6.405 1.98e-10 ***
## dew_point_temperature_c  3.819e-02  9.051e-03   4.219 2.59e-05 ***
## solar_radiation_mj_m2   -4.150e-03  3.073e-02  -0.135 0.892598    
## rainfall_mm             -3.069e-01  1.014e-01  -3.025 0.002527 ** 
## snowfall_cm             -7.377e-02  2.196e-02  -3.360 0.000799 ***
## holiday1                -4.533e-01  3.432e-02 -13.210  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 24.48278)
## 
##     Null deviance: 147711  on 1627  degrees of freedom
## Residual deviance:  41915  on 1595  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
y_pre_train <- predict(md.qps.wt,trainset,type = "response")
y_pre_test <- predict(md.qps.wt,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 83.3299356536891"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 95.7229074446783"
myGLM = function(formula, data) {
  glm(formula, data, family = quasipoisson(link = log))
}

myPredictGLM = function(object, newdata){
  predict(object, newdata , type="response")
}

md.qps.cv.wt <- errorest(rented_bike_count ~ .,
                         data = data.model.wt, predict  = myPredictGLM,model = myGLM,estimator="cv",
                     est.para=control.errorest(k=5, predictions = TRUE))
md.qps.cv.wt
## 
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model.wt, 
##     model = myGLM, predict = myPredictGLM, estimator = "cv", 
##     est.para = control.errorest(k = 5, predictions = TRUE))
## 
##   5-fold cross-validation estimator of root mean squared error
## 
## Root mean squared error:  87.8302

Summer

data.model.sm <- data.model[ data.model$ seasons =="Summer",- which(colnames(data.model) %in%  c('seasons'))]
rs <- slit_test_train(data.model.sm,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])


md.qps.sm <- glm(rented_bike_count ~  .,family = quasipoisson('log'),data=trainset )
summary(md.qps.sm)
## 
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson("log"), 
##     data = trainset)
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              9.574e+00  2.559e-01  37.416  < 2e-16 ***
## hour1                   -3.120e-01  4.996e-02  -6.246 5.41e-10 ***
## hour2                   -6.352e-01  5.509e-02 -11.530  < 2e-16 ***
## hour3                   -1.063e+00  6.423e-02 -16.548  < 2e-16 ***
## hour4                   -1.468e+00  7.467e-02 -19.660  < 2e-16 ***
## hour5                   -1.357e+00  7.392e-02 -18.361  < 2e-16 ***
## hour6                   -7.258e-01  5.821e-02 -12.467  < 2e-16 ***
## hour7                   -3.902e-02  4.723e-02  -0.826   0.4088    
## hour8                    2.510e-01  4.454e-02   5.636 2.06e-08 ***
## hour9                   -1.987e-01  5.042e-02  -3.941 8.48e-05 ***
## hour10                  -4.527e-01  5.529e-02  -8.187 5.43e-16 ***
## hour11                  -3.784e-01  5.618e-02  -6.735 2.28e-11 ***
## hour12                  -2.834e-01  5.707e-02  -4.966 7.59e-07 ***
## hour13                  -2.695e-01  5.909e-02  -4.561 5.49e-06 ***
## hour14                  -2.321e-01  5.683e-02  -4.085 4.64e-05 ***
## hour15                  -8.945e-02  5.381e-02  -1.662   0.0967 .  
## hour16                   9.731e-02  5.075e-02   1.917   0.0554 .  
## hour17                   4.401e-01  4.573e-02   9.623  < 2e-16 ***
## hour18                   8.409e-01  4.106e-02  20.479  < 2e-16 ***
## hour19                   7.612e-01  4.097e-02  18.580  < 2e-16 ***
## hour20                   7.137e-01  4.045e-02  17.647  < 2e-16 ***
## hour21                   6.627e-01  4.043e-02  16.391  < 2e-16 ***
## hour22                   5.279e-01  4.053e-02  13.026  < 2e-16 ***
## hour23                   2.318e-01  4.323e-02   5.361 9.47e-08 ***
## temperature_c           -8.505e-02  9.592e-03  -8.867  < 2e-16 ***
## humidity_percent        -2.089e-02  2.856e-03  -7.314 4.09e-13 ***
## wind_speed_m_s           5.763e-03  8.675e-03   0.664   0.5066    
## visibility_10m          -8.665e-05  1.548e-05  -5.598 2.54e-08 ***
## dew_point_temperature_c  5.590e-02  1.015e-02   5.509 4.20e-08 ***
## solar_radiation_mj_m2    9.602e-02  1.674e-02   5.735 1.16e-08 ***
## rainfall_mm             -1.075e+00  8.136e-02 -13.217  < 2e-16 ***
## snowfall_cm                     NA         NA      NA       NA    
## holiday1                -3.540e-02  4.355e-02  -0.813   0.4165    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 68.63507)
## 
##     Null deviance: 653959  on 1617  degrees of freedom
## Residual deviance: 116143  on 1586  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
y_pre_train <- predict(md.qps.sm,trainset,type = "response")
y_pre_test <- predict(md.qps.sm,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 287.701859290571"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 320.851827259716"
myGLM = function(formula, data) {
  glm(formula, data, family = quasipoisson(link = log))
}

myPredictGLM = function(object, newdata){
  predict(object, newdata , type="response")
}

md.qps.cv.sm <- errorest(rented_bike_count ~ .,
                         data = data.model.sm, predict  = myPredictGLM,model = myGLM,estimator="cv",
                     est.para=control.errorest(k=5, predictions = TRUE))
md.qps.cv.sm
## 
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model.sm, 
##     model = myGLM, predict = myPredictGLM, estimator = "cv", 
##     est.para = control.errorest(k = 5, predictions = TRUE))
## 
##   5-fold cross-validation estimator of root mean squared error
## 
## Root mean squared error:  298.6332

Sau khi phân chia ra các mùa cụ thể ta thấy dữ liệu trong từng nhóm đồng chất hơn nên độ phân tán mô hình thấp hơn và cũng như rmse cũng thấp mô hình tổng.

Splines

Cải thiện cho biến humidity_percent và temperature_c là 2 biến tác động lớn trong mô hình hồi quy poisson

rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])
knots_h <- quantile(trainset$humidity_percent, probs = c(0.25, 0.5, 0.75)) 
knots_t <- quantile(trainset$temperature_c, probs = c(0.25, 0.5, 0.75)) 
md.ps.sl <- glm(rented_bike_count ~  hour + bs(temperature_c,knots_t,3) + bs(humidity_percent,knots_h,3) + wind_speed_m_s + dew_point_temperature_c +  solar_radiation_mj_m2 + rainfall_mm + snowfall_cm + holiday + seasons,
                           data = trainset,family = quasipoisson)
y_pre_train <- predict(md.ps.sl,trainset,type = "response")
y_pre_test <- predict(md.ps.sl,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 253.687888119918"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 259.658732815563"
print(paste("Độ phân tán:",sum(residuals(md.ps.sl,"pearson" ) ^2 ) /md.ps.sl$df.residual ))
## [1] "Độ phân tán: 73.9049528539409"

Mô hình có sự tương tác

rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])


md.ps.int <- glm(rented_bike_count ~ (hour + temperature_c + humidity_percent + wind_speed_m_s + dew_point_temperature_c +  solar_radiation_mj_m2 + rainfall_mm + snowfall_cm + holiday) * seasons,
                           data = trainset,family = quasipoisson)
y_pre_train <- predict(md.ps.int,trainset,type = "response")
y_pre_test <- predict(md.ps.int,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))
## [1] "Rmse của trainset: 247.393906091801"
print(paste("Rmse của testset:",rmse_test))
## [1] "Rmse của testset: 253.290239415538"
print(paste("Độ phân tán:",sum(residuals(md.ps.int,"pearson" ) ^2 ) /md.ps.int$df.residual ))
## [1] "Độ phân tán: 70.4196996642662"
anova(md.ps.int, test="Chisq")
## Analysis of Deviance Table
## 
## Model: quasipoisson, link: log
## 
## Response: rented_bike_count
## 
## Terms added sequentially (first to last)
## 
## 
##                                 Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                             6296    3247619              
## hour                            23  1248049      6273    1999569 < 2.2e-16 ***
## temperature_c                    1   963614      6272    1035955 < 2.2e-16 ***
## humidity_percent                 1    46470      6271     989485 < 2.2e-16 ***
## wind_speed_m_s                   1     2585      6270     986900 1.380e-09 ***
## dew_point_temperature_c          1     4947      6269     981953 < 2.2e-16 ***
## solar_radiation_mj_m2            1     1586      6268     980367 2.076e-06 ***
## rainfall_mm                      1    72584      6267     907783 < 2.2e-16 ***
## snowfall_cm                      1    12329      6266     895454 < 2.2e-16 ***
## holiday                          1     9547      6265     885907 < 2.2e-16 ***
## seasons                          3   245323      6262     640584 < 2.2e-16 ***
## hour:seasons                    69    75782      6193     564802 < 2.2e-16 ***
## temperature_c:seasons            3   122509      6190     442293 < 2.2e-16 ***
## humidity_percent:seasons         3     3311      6187     438982 3.468e-10 ***
## wind_speed_m_s:seasons           3     1568      6184     437414 5.748e-05 ***
## dew_point_temperature_c:seasons  3     1461      6181     435953 0.0001194 ***
## solar_radiation_mj_m2:seasons    3      480      6178     435473 0.0778307 .  
## rainfall_mm:seasons              3      470      6175     435003 0.0828562 .  
## snowfall_cm:seasons              1      936      6174     434066 0.0002666 ***
## holiday:seasons                  3     3235      6171     430831 5.878e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mô hình khi có sự tương tác giữa seasons với các yếu tố còn lại thực sự hiểu quả hơn, tuy nhiên độ phức tạp cũng tăng.

Co hệ số

X <- data.numeric |> select(-rented_bike_count) |> scale() |> as.matrix()
y <- data.numeric |> select(rented_bike_count)  |> unlist() 
out.cv.lasso <- cv.glmnet(x = X, y = y, alpha = 1,
type.measure = "mse", nfolds = 5,
family = poisson)
print(out.cv.lasso)
## 
## Call:  cv.glmnet(x = X, y = y, type.measure = "mse", nfolds = 5, alpha = 1,      family = poisson) 
## 
## Measure: Mean-Squared Error 
## 
##     Lambda Index Measure   SE Nonzero
## min  0.421    74   99637 2447      35
## 1se  6.257    45  101888 2333      29
beta.lambda.lasso <- out.cv.lasso$lambda.min
md.lasso <- glmnet(x = X, y = y, alpha = 1, family = quasipoisson)
predict(md.lasso, s = beta.lambda.lasso, type = "coefficients")
## 39 x 1 sparse Matrix of class "dgCMatrix"
##                                    s1
## (Intercept)              6.261395e+00
## hour.0                  -1.526357e-02
## hour.1                  -6.417294e-02
## hour.2                  -1.293567e-01
## hour.3                  -2.001664e-01
## hour.4                  -2.892629e-01
## hour.5                  -2.706017e-01
## hour.6                  -1.218734e-01
## hour.7                   2.201748e-02
## hour.8                   1.125281e-01
## hour.9                   7.805370e-03
## hour.10                 -4.903058e-02
## hour.11                 -3.680459e-02
## hour.12                 -1.586947e-02
## hour.13                 -1.882491e-02
## hour.14                 -1.817354e-02
## hour.15                  .           
## hour.16                  3.154337e-02
## hour.17                  8.232566e-02
## hour.18                  1.489885e-01
## hour.19                  1.145066e-01
## hour.20                  1.006405e-01
## hour.21                  1.008192e-01
## hour.22                  8.100793e-02
## hour.23                  2.202673e-02
## temperature_c            1.116877e-01
## humidity_percent        -1.881183e-01
## wind_speed_m_s          -7.145285e-03
## visibility_10m          -5.719105e-03
## dew_point_temperature_c  2.401607e-01
## solar_radiation_mj_m2    2.427883e-02
## rainfall_mm             -3.306467e-01
## snowfall_cm             -1.845556e-02
## seasons.Autumn           8.104036e-02
## seasons.Spring           .           
## seasons.Summer           .           
## seasons.Winter          -3.564635e-01
## holiday.0                3.714379e-02
## holiday.1               -6.736188e-10

Với phương pháp co hệ số ta co thể loại bỏ một số biến không tác động đến mô hình.

Tổng kết:

Sau thực hiện mô hình ta rút ra được những insight, đó là sự biến đổi số lượng xe thuê giá phụ thuộc vào mùa, thời gian,và sự tác động các biến giải thích đến số lượng xe thuê. Từ đấy, phục vụ cho các mục đích kinh doanh, marketing,…